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Continuous-time random walk (CTRW) is a model of anomalous sub-diffusion in which particles 
are immobilized for random times between successive jumps. A power-law distribution of the waiting 



times, ip(r) 



-(l+oO 



, leads to sub-diffusion ((s 2 ) ~ t a ) for < ct < 1. In closed systems, the 



long stagnation periods cause time-averages to divert from the corresponding ensemble averages, 
which is a manifestation of weak ergodicity breaking. The time-average of a general observable 
U = j f* U[x(t)]cIt is a functional of the path and is described by the well known Feynman-Kac 
equation if the motion is Brownian. Here, we derive forward and backward fractional Feynman-Kac 
equations for functionals of CTRW in a binding potential. We use our equations to study two 
specific time-averages: the fraction of time spent by a particle in half box, and the time-average of 
the particle's position in a harmonic field. In both cases, we obtain the probability density function 
of the time-averages for t — >• oo and the first two moments. Our results show that indeed, both the 
occupation fraction and the time-averaged position are random variables even for long-times, except 
for q = 1 when they are identical to their ensemble averages. Using the fractional Feynman-Kac 
equation, we also study the dynamics leading to weak ergodicity breaking, namely the convergence 
of the fluctuations to their asymptotic values. 



PACS numbers: 05.40.Fb,05.40.Jc,05.10.Gg,02.50.Ey 



I. INTRODUCTION 

The time-average of an observable U(x) of a diffusing 
particle is defined as 



u= 1 - 



U[x{T)]dT, 



(1) 



where x(t) is the particle's trajectory. For Brownian mo- 
tion in a binding potential V(x) and in contact with a 
heat bath, ergodicity leads to 



U(x)G eq (x)dx, 



(2) 



where G cq {x) = e ~v(x)l{k B T) j Z ig Boltzmann distribu- 
tion and (E/)th * s the thermal average. The equality of 
time- and ensemble averages in ergodic systems is one of 
the basic presuppositions of statistical mechanics. 

In the last decades it was found that in many systems, 
the diffusion of particles is anomalously slow, as charac- 
terized by the relation (a; 2 ) ~ t a with < a < 1 0- 
Anomalous sub-diffusion is commonly modeled as a 
continuous-time random walk (CTRW): nearest- neighbor 
hopping on a lattice, with waiting times betweeniumps 
distributed as a power-law with infinite mean [H, |(3|. 

For closed systems, the long immobilization periods of 
CTRW result in deviation of time-averages from ensem- 
ble averages even for long times @-[l(j. Although there 
are no inaccessible regions in the phase space (i.e., there 
is no strong ergodicity breaking), the divergence of the 
mean waiting time results in some waiting times of the 
order of magnitude of the entire experiment. Therefore, 
a particle does not sample the phase space uniformly in 
any single experiment, leading to weak ergodicity break- 
ing 0. 



Two examples of particularly interesting time- 
averages, which we study in this paper, are given be- 
low. For a particle in a bounded region, the occupa- 
tion fraction is defined as A = \ Q[x(r)]dT, namely, 
it is the fraction of time spent by the particle in the 
positive side of the region 

[Hill- Generally, the occu- 
pation fraction can be defined for any given subspace. 
Consider, for example, a particle in a sample illumi- 
nated by a laser, where the particle emits photons only 
when it is under the laser's focus. The occupation frac- 
tion is proportional to the total emitted light [H], EH- 
Next, the time-averaged position of a particle is defined 
as x — j L x(r)dT. Recent advances in single particle 
tracking technologies enable the experimental determina- 
tion of the time-average of the position of beads in poly- 
mer networks [TfJ [ItJ and of biological macro-molecules 
and small organelles in living cells |T8l - (2l| . Since in many 
physical and biological systems the diffusion is anoma- 
lous, the study of occupation fractions or time-averaged 
positions in sub-diffusive processes such as CTRW is of 
current interest. 

Time-averages are closely related to functionals, which 
are defined as A = J U[x(r)\dT and have many applica- 
tions in physics, mathematics and other fields [22j |. De- 
note by G(x, A, t) the joint PDF of finding, at time t, the 
particle at x and the functional at A. The Feynman-Kac 
equation states that for a free Brownian particle [23[ : 



d_ 

at 



G(x, P ,t) 



d 2 

K l-g^2 G (x,p,t) 



■pU{x)G(x,p,t), (3) 



where G(x,p,t) is the Laplace transform A — > p of 
G(x, A, t) and K\ is the diffusion coefficient. Recently, we 
developed a fractional Feynman-Kac equation for anoma- 
lous diffusion of free particles [HI, [25j] . As time-averages 
are in fact scaled functionals: U = A/t, a generalized 



2 



Feynman-Kac equation for anomalous functionals in a 
binding field would be invaluable for the study of weak 
crgodicity breaking. Currently, no such equation exists 
and weak ergodicity breaking was investigated only in the 
t — > oo limit or using functional- and potential-specific 
methods 0-Q3- 

In this paper, we obtain an equation for functionals of 
anomalous diffusion in a force field F(x). The equation 
takes the following form (reported without derivation in 
0): 



d_ 

at 

K, 



G(x,p,t) = 

d 2 d F{x) 
dx 2 dx ksT 



(4) 

V 1 t - a G(x,p,t)- P U(x)G(x,p,t). 



The symbol D* - Q is a fractional substantial derivative, 
equal in Laplace t — > s space to [s + pU(x)] 1 ~ a p6l. [27j. 
and K a is a generalized diffusion coefficient. Solving Eq. 
(H|) for G(x,p,t), inverting p — > A and integrating over 
all x yields G(A, t), the PDF of A at time t. Changing 
variables A —> A/t = U, one finally comes by G(U,t), 
the (time-dependent) PDF of U. Weak crgodicity break- 
ing can then be determined by looking at the long-times 
properties of G(U, t): iiU is not identically equal to (U) th 
for t — > oo, ergodicity is broken. Moreover, if G(U,t) or 
the moments of U can be found also for t < oo, the ki- 
netics of weak ergodicity breaking can be uncovered. 

In the rest of the paper, we derive Eq. (j4j) as well as a 
backward equation and an equation for time-dependent 
forces. We then apply our equation to the two examples 
given above: the occupation fraction in a box and the 
time-averaged position in a harmonic potential. In both 
cases, we calculate the long-times limit of G(U, t) and the 

fluctuations ((AC/) 2 ) = (jj 2 ^ - (JJ) 2 . Wc demonstrate 

that for sub-diffusion both systems exhibit weak ergod- 
icity breaking, and that the fluctuations decay as t~~ a to 
their asymptotic limit. Part of the results for the fluctua- 
tions of the time-averaged position were briefly reported 
in M- 



II. DERIVATION OF THE FRACTIONAL 
EQUATIONS 

A. The forward equation 

1. Continuous-time random walk 

In the continuous-time random walk model, a particle 
is placed on an one-dimensional lattice with spacing a 
and is allowed to jump to its nearest neighbors only. The 
probabilities of jumping left L(x) and right R(x) depend 
on F(x), the force at x (see next subsection for derivation 
of these probabilities). If F(x) = 0, then R(x) — L(x) = 
1/2. Waiting times between jump events are independent 
identically distributed random variables with PDF ip(r), 



and are independent of the external force. The initial 
position of the particle, xq, is distributed according to 
Gq(x). The particle waits in xq for time r drawn from 
i/}(t), and then jumps to either xq + a (with probability 
R{x)) or xq — a (with probability L(x)), after which the 
process is renewed. We assume that the waiting time 
PDF scales as 



B a 



|r(-a)| 



-(1+a) 



(5) 



where < a < 1. With this PDF, the mean waiting time 
is infinite and the process is sub-diffusive: for F(x) = 
0, xp = 0, and for an infinite open system, (a: 2 ) ~ t a 
[281 ] . We also consider the case when the mean waiting 
time is finite, e.g., an exponential distribution ip(r) = 
e — T /( r ) j ( r )_ This leads to normal diffusion (x 2 ) ~ t and 
we therefore refer to this case as a = 1. For discussion 
on the effect of an exponential cutoff on Eq. , see [2!| . 
Below, we derive the differential equation that describes 
the distribution of functionals in the continuum limit of 
this model. 



2. Derivation of the equation 

Define A = ^U[x{r)]dr and define G(x, A, t) as the 
joint PDF of x and A at time t. For the particle to be at 
(x, A) at time t, it must have been at [x,A — tU(x)] at 
time t — t when the last jump was made. Let x( x i A, t)dt 
be the probability of the particle to jump into (x, A) in 
the time interval [t, t + dt]. We have, 



G{x,A,t) 



W{r)x[x,A-TU(x),t-T]dT, (6) 



where W(t) = 1 — i/ , ( r ')^ r ' i s the probability for not 
moving in a time interval of length r. 

To calculate \i note that to arrive to (x, A) at time 
t, the particle must have arrived to either [x — a, A — 
tU (x — a)] or [x + a, A — tU (x + a)] at time t — r when 
the previous jump was made. Therefore, 

X (x,A,t) = G (x)S(A)5(t) (7) 

+ / 4>(t)L(x + a)x[x + a, A - tU(x + a),t - r]dr 



+ / ^P(t)R(x — a)x[x — a, A — tII(x — a), t — r]dr. 
Jq 

The term Go(x)S(A)S(t) corresponds to the initial con- 
dition, namely that at t = 0, A = and the particle's 
position is distributed as Gq(x). 

Assume that U(x) > for all x and thus A > 
(an assumption we will relax in Section III A 3[) . Let 
X{x,p,t) = f °° e~ pA x(x,A,t)dA be the Laplace trans- 
form A — > p of x(x, .A, t) (we use along this work the 
convention that the variables in parenthesis define the 
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space we are working in). Laplace transforming Eq. ([?]) 
from A to p, we find 



X (x,p,t) = G (x)6(t) 

+ L(x + a) [ i>(T)e- pTUix+a) x(x + a,p,t-T)d7 



(8) 



+ R(x-a) / ^(r)e- prC/ ^- a )x(a;-a,p,i-r)rfr. 

Laplace transforming Eq. (|5J| from i to s using the con- 
volution theorem, 

X (x,p,s)=Go(x) (9) 
+ L(x + a)ip[s + pU (x + a)]x(x + a,p, s) 
+ R(x — a)tp[s + pU (x — a)]x(x — a,p, s), 

where ip(s) is the Laplace transform of the waiting time 
PDF. Let x(k, p, s) — e lkx x(x, p, s)dx be the Fourier 
transform x — > k of x- Fourier transforming Eq. ([9]) and 
changing variables x ± a — > x, 



X{k,p,s) = G (k) 



(10) 



-ika 



ika 



5 L(x)^[s + pU(x)]x(x,p, s)dx 



e lkx R(x)iJj[s +pU(x)]x(x 1 p,s)dx, 



where Go(fc) is the Fourier transform of the initial con- 
dition. 

We now express L{x) and R(x) in terms of the po- 
tential V(x). Assuming the system is coupled to a heat 
bath at temperature T and assuming detailed balance, 
we have fiol |28| 



L(x) exp 



V(x) 



knT 



R(x — a) exp 



V(x - a) 



knT 



(11) 



If the lattice spacing a is small we can expand 

aF{x) 



exp 



V(x-a)~ 




r v( x y 




k B T 


« exp 


k B T _ 





1 



knT 



+ 0(a 2 ) 
(12) 



where we used F(x) = — V'(x). Expanding R(x) and 
L(x) for « 
1/2 for F(x) = 0, 



L(x) for r§^jl <C 1, using the fact that R(x) — L(x) 



R(x) 



l + c 



aF{x) 



knT 



1-L(x), 



(13) 



where c is a constant to be determined. Combining Eqs. 
([TTj). (|T2"j) , and ([TB")) , we have, up to first order in a 



1 - c 



aF(x) 
k B T 



[ a^(x)" 






fc B T _ 







This gives, again up to first order in a, c = 1/2. We can 
thus write, 



l 



aF(x) 
2k B T 



L(x) « - 
V y 2 



2k B T 



(14) 



Substituting Eq. dT4J) in Eq. (|10|). we obtain, 

X(fc,p, s) w G (fc) 
1 



2 6 
1 



-ika 



2 
1 

2 6 ~ 



e lfc3: ^[s +pC/(x)]x(a;,p, s)cte 



e ifcc ^r^V' [s + PU (x)] X (x, p, s)dx 
2k B T 

ijj[s + pU(x)]x(x,p, s)dx 



, ikx 



Aha 



e lkx ^^ip[s+pU(x)}x(x,p,s)dx. 
2k B l 



Applying the Fourier transform identity iF{x/(x)} 

l dk- 



-i-Srf{k), the last equation simplifies to 



X (k,p, s)«G (fc)+ 



a F (-id 
cos(fca) + i sin(fca) — 



2fc B T 
X(k,p,s). (15) 



The symbols F (— and U (—i-j^) represent the orig- 
inal functions F{x) and U(x), but with — i-^g as their 
arguments. Note that the order of the terms is im- 
portant: for example, cos(/ca) does not commute with 
ip [s + pU (— i-§^)\ ■ The formal solution of Eq. (fi~5|) is 



1 - 



cos(fca) + i sin(fca)- 



aF -i 



t)k > 



2k B T 







s + pU —i 



dk 



Go(fc). 



(16) 



We next use our expression for x to calculate G(x, A, t). 
Transforming Eq. (JB|) (x, A, t) — > (k,p, s), 

^p,^ 1 -^;^.^ ^^), d7) 

s +p u hi) 

where we used the fact that W(s) = [1 — ip(s)]/s. Sub- 
stituting Eq. (|16[) into (fTT|) . we have 



G(fc,p, s) 



1 - 



(18) 



cos(fca) + zsin(fca) 



2k B T 



s + pU —i 



dk 



G (k). 



To derive a differential equation for G(x,p, t), we use 
the small s expansion of tp{s). For < a < 1, where the 
waiting time PDF is ip{r) ~ 3^-^+^ /W(-a)\ (Eq. 
((SJ)), the Laplace transform for small s is [3] 

V>(s) ~ 1 — B a s a . (19) 
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The case a = 1 is also described by Eq. ([T9|) . if we 
identify B\ with the mean waiting time (t) . Substituting 
Eq. ((Ill) in Eq. fl8]) . and using cos(A:a) w 1 - fc 2 a 2 /2 
and sin(fca) « fca, we obtain 



G(k,p,s) 



s+pU 



dk 



(20) 



ik- 



k B T 



This is simply the (integer) Feynman-Kac equation (|3|), 
extended to a general force field F(x). 

The fractional Fokker-Planck equation. — For p = 0, 
G(x,p = Q,t) = J °° G(x,A,t)dA reduces to G(x,t), 
the marginal PDF of finding the particle at x at time 
t regardless of the value of A. Correspondingly, Eq. 
(l22l) reduces to the fractional Fokker-Planck equation 



G (k), 



where we defined the generalized diffusion coefficient [2£ 

,,2 



K a — lim 



a 2 ,B Q -s.O 2B n ' 



(21) 



with units m 2 /scc Q . Rearranging the expression in Eq. 

(2D), 

sG(k,p, s) - G (k) = - P U G{k,p, s) 



k 2 - ik~^ — ^ 



k B T 



s + pU —i 



dk 



l-a 



G(k,p, s). 



Inverting k — > x, s t, we finally obtain our fractional 
Feynman-Kac equation: 

9 



dt 



G(x,p,t) = K a L FP T)l- a G(x,p,t) - P U(x)G(x,p,t). 

(22) 

The symbol Lfp represents the Fokker-Planck operator, 
d 2 d F(x) 



FP 



dx 2 dx k B T ' 



(23) 



and the initial condition is G(x, A, t = 0) = Gq(x)S(A), 
or G(x,p,t = 0) = Go(x). The symbol D]~ a represents 
the fractional substantial derivative operator introduced 
in [Hill: 

£{'Dl- a G(x,p,t)} = [ S +pU(x)] 1 - a G(x,p,s), (24) 

where L{f(t)} = J °° e~ st f(t)dt is the Laplace transform 
t — > s. In t space, 



Vl- a G(x,p,t) 



(25) 



T(a) 



d_ 

dt 



P U(x) 



* p -(t-r)pU(x) 



— r-j G(x, p, T)dT. 



Thus, due to the long waiting times, the evolution of 
G{x,p,t) is non-Markovian and depends on the entire 
history. 



d_ 



G{x,t) = K a Z FV 'Dl- L a t G{x,t), 



(27) 



^ere D 1 ^ 



■ Several 



Dl a \ „ is the Riemann-Liouville 

1 lp=0 

fractional derivative operator. In Laplace s space, 
V 1 K £G(x,s) = s 1 - a G(x,s). 

Free particle. — For F(x) = 0, £pp 
applications of this special case were treated in [2 

A general functional. — When the functional is not 
necessarily positive, the Laplace transform A — > 
p is replaced by a Fourier transform G(x,p,t) = 
j^' ao e lpA G{x,A,t)dA. The fractional Feynman-Kac 
equation looks like (|22|) . but with p replaced by — ip, 







—G(x,p,t) 



if Q £ FP D t i - Q G(x,p, t) + i P U{x)G(x,p, t), 

(28) 

where D] a — > [s — ^^(a;)] 1 a in Laplace s space. The 
derivation of Eq. (j28j) is similar to that of (|22|) (see 
for more details). 

Time- dependent force. — Anomalous diffusion with a 
time-dependent force is of recent interest [33l - l37j . When 
the force is time-dependent, we assume the probabilities 
of jumping left and right are determined by the force 
at the end of the waiting period [H, As we show 
in Appendix A, the equation for the PDF G(x,p,t) is 
similar to Eq. (|2"2"j) : 

^G{x,p,t) = K a L^lv 1 t - a G(x lP ,t)-pU(x)G(x,p,t), 

(29) 

but where 



(t) = 2 d F(x,t) 
FP 8x 2 dx k R T 



is the time-dependent Fokker-Planck operator. For p = 
0, Eq. ([29| reduces to the recently derived equation for 
the PDF of x HI. 



B. The backward equation 



3. Special cases and extensions 

Normal diffusion. — For a = 1, or normal diffusion, 
the fractional substantial derivative equals unity and we 
have 

^-G(x,p,t)=K 1 L FF G(x,p,t)-pU{x)G{x,p,t). (26) 



The forward equation describes G(x,A, t), the joint 
PDF of x and A. Consequently, if we are interested only 
in the distribution of A, we must integrate G over all x, 
which could be inconvenient. We therefore develop below 
an equation for G xo (A, t) — the PDF of A at time t, given 
that the process has started at x . This equation, which 
is called the backward equation, turns out very useful 
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in practical applications (see, e.g., [22|, [25| and Section 

ttYU). 

According to the CTRW model, the particle starts at 
x = .To and jumps at time t to either xq + a or xq — a. 
Alternatively, the particle does not move at all during 
the measurement time [0,i]. Hence, 

G Xo (A,t) = W(t)5[A-tU(x )] (30) 

+ / 1p(T)R(x )G Xo+a [A-TU(x ),t-T)dT 

Jo 

+ / i){T)L{x Q )G Xo - a {A-TU(xo),t-T\dr. 
Jo 

Here, tU(xq) is the contribution to A from the paus- 
ing time at xq in the time interval [0,r]. The first term 
on the rhs of Eq. (|3T))) describes a motionless particle, 
for which Alt) = tlf(xo)- We now transform Eq. (|30)) 
(xo,A, t) — > (ko,p,s), using techniques similar to those 
used in Section III A 21 In the continuum limit, a — > 0, 
this leads to, 



Gk (p, s) 



1-$ 




s+pU 







S(k 



cos(fcoa) 



S+pU[ ~ l dk- 



2k B T 



i sin(fcoa) 



Gk {p,s). 



We then expand ip(s) « 1 — B a s a , cos(fcoa) ~ 1 — k 2 a 2 /2, 
and sin(fc a) ~ k$a. After some rearrangements, 

sG ko (p,s) - 6(k ) = -pU ( -i-^r- ] G ko (p,s) 



Ka 



s + pU —i 



d_ 

dkr 



feo H ; — - , * fc o 



<9fc 



G ko {v,s). 



Inverting fco — > xo and s — > i, we obtain the backward 
fractional Feynman-Kac equation: 



d_ 

dt' 

where 



G Xo (p,t) = K a ^ 1 t - a L^G X0 { P ,t)-pU(x )G X0 (p,t), 

(31) 

(32) 



(B) _ 
FP 



d 2 F(x ) d 



Oxq ksT dxo 



is the backward Fokkcr-Planck operator. The initial con- 
dition is G X0 (A,t = 0) = 8(A), or G xo (p,t = 0) = 1. 

Note the (+) sign of £p P and the order of the operators 
in its second term, which are opposite to those of £fp 



(Eq. (|23|0 . Here, T)\~ a equals in Laplace t s space 
[s + pU(xq)} 1 ^ 01 . In Eq. (|2"2")l the operators depend on x 
while in Eq. (|3T|) they depend on xo . Therefore, Eq. (|22|) 
is a forward equation while Eq. pip is a backward equa- 
tion. Notice also that in Eq. (|3"Tj) , the fractional deriva- 
tive operator appears to the left of the Fokker-Planck 
operator, in contrast to the forward equation (122ft . 



III. THE PDF OF U FOR LONG TIMES 

For long measurement times, it is possible to use the 
fractional Feynman-Kac equation to obtain an expression 
for the PDF of a general time-average: 



A 

7' 



fiU[x(T)]dT 
t 



We write first the forward equation ([22]) in Laplace s 
space: 



[s + pU(x)\ G(x,p, s) — G (:c 
d 2 d F(x) 



(33) 



dx 2 dx k B T 



[s+pU(x)Y- a G{x,p lS ). 



CTRW functionals scale linearly with the time, A ^ t, 
and therefore, as shown in [38|, G(p, s) = g(p/s)/s, where 
g is a scaling function. Since we are interested in the 
t — > oo limit, we take s and p to be small, with their ratio 
finite. We therefore expect G(x,p,s) ~ s _1 (indeed, see 
Eq. (|36[) below), and consequently, both terms on the 
lhs of (|3"3")l scale as s°. However, the rhs of (j3"3"|) scales as 
s~ a , and therefore for small s the lhs is negligible. The 
forward equation thus reduces to 



K a 



c) 2 



d F(x) 



s+pU{x)] l - a G(x,p,s) =0. 



dx 2 dx k B T 
The solution of the last equation is 

G(x,p,s) = C(p, s)[s + pUix)]"- 1 exp 



k R T 



(34) 



where G(p, s) is independent of x. To find G, we integrate 
Eq. ([33j) over all x: 



[s + p£/(x)] G(x,p, s)cfe -1 = 0, 



(35) 



which is true, because for a binding field, G(x,p,s) and 
its derivative vanish for large \x\. Substituting G from 
Eq. {35 into Eq. {35} gives 



C(p,s) = \[ [ S +pU(x)] a e X p 
Therefore, 



fc R T 



f/.r 



G(x,p, s) 



[s+pC/(a:)] Q " 1 exp 


k B T 




exp 


V(x) 
k B T _ 


dx 



(36) 
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Integrating Eq. (|36[) over all x, 



The distribution 



G(p,s) 



cxp 



YM 

k B T 



dx 



YM 

k B T 



dx 



(37) 



where G(p, s) is the double Laplace transform of G(A, t), 
the PDF of A at time t. The last equation is the con- 
tinuous version of the result derived using a different ap- 
proach in [j| [To| ■ As in_p, [Toj , Eq. ([37| can be inverted, 
using the method of [38[, to give the PDF of U = A/t, 



G(U) = Sin(?m) x 



(38) 



[I^iU)] 2 + [I<(UW+2co S (7r a )I>(U)I<(Uy 



where 



cxp 



U<U(x) 



and 



I*(U)= l_ exp 

IU>U(x) 



YM 



YM 



[U{x)-U\ a dx 



\U-U{x)] a dx. 



For normal diffusion, a = L, the PDF is a delta func- 
tion G(U) = 5 [U- (U) th \ For anomalous sub- 
diffusion, a < 1, U is a random variable, different from 
the ensemble average. This behavior of the time-average 
results from the weak ergodicity breaking of the sub- 
diffusing system. Similar results hold when U(x) is not 
necessarily positive: the Laplace transform A — >• p is 
replaced by a Fourier transform and in Eq. (|37j) . p is 
replaced by — ip. 



IV. APPLICATIONS: WEAK ERGODICITY 
BREAKING 

In this section we present two applications of the frac- 
tional Feynman-Kac equation: the occupation fraction 
in a box and the time-averaged position in a harmonic 
potential. We demonstrate weak ergodicity breaking in 
both cases and investigate the convergence to the asymp- 
totic limits. 



A. The occupation fraction in the positive half of a 

box 



We study the problem of the occupation time in x > 
for a sub-diffusing particle moving freely in the box 



extending between [- 



L L 
2 ' 2 . 



0,1 EH. 



Define the occupation time in x > as T + = 
J*Q[x{T)]dT (namely U{x) = Q(x)). To find the PDF 
of T + , we write the backward fractional Feynman-Kac 
equation (|3"Tj) in Laplace s space: 



sG X0 (p, s) - 1 



(39) 
x < 0, 

The equation (|39p is subject to the boundary conditions: 



'K a s l -<*-§^G Xo {p,s) 
K a (s+p) 1 - a £ J G X0 (p,s)~pG X0 { P ,s) x > 0. 



UX x =± . 

The solution of the last equation is 
Cq cosh 



= 0. 



Gx (p,s) 



C\ cosh 



(f+p) 



s+p 



x < 0, 
x > 0. 



(40) 

Matching G and its derivative at xq = gives the equa- 
tions: 



Cq cosh 



Ls a l 2 



1 



Ci cosh 



L(s + p) a /' 2 



2^K a 



1 



C s a/2 sinh ( i^ll ) = -d(s + p) Q/2 sinh 



2^ 



s+p 
L{s+p) a / 2 ~ 
2JKZ 



Solving these equations for Co and C\ and substituting 
x Q = in Eq. (|4T)j) gives, after some algebra, 

G (p,s)= (41) 

W 2 - 1 tanh [(sr) Q / 2 ] + (s + p)"/ 2 " 1 tanh {[r(s + p)]" /2 } 

W 2 tanh [(sr) Q / 2 ] + (s + p) Q / 2 tanh { [r(s + p)] Q/2 } 

where we defined r Q = L 2 /(4:K a ). This equation was 
previously derived in Q using a different method. Eq. 
(j4Tj) describes the PDF of T+ for all times, but cannot 
be directly inverted. For long times, or (sr) Q / 2 <C 1, 



G (p,s) 



+ (s+p) 



a-l 



s a + (s+p) a 



(42) 



This can be inverted to give the PDF of A = T + /t, or the 
occupation fraction [H, [38| , 



G(A) 



sin(7ra) 



A" _1 (l - A) c 



7T A 2a + (l-A) 2Q + 2cos(7ra)A a (l-A) Q ' 

(43) 

Eq. 03} is called Lamperti's PDF |39j. Note that Eqs. 
(|42|) and (|43|) can also be derived directly from the gen- 
eral long-times limit, Eqs. (|3"T|) and (J3HJ) , respectively. 
Whereas the PDF of the occupation fraction for a free 
particle is also Lamperti's [25| , in the free particle case 
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FIG. 1: The PDF of the occupation fraction in half-space 
for a particle in the box [—1, 1]. CTRW trajectories were 
generated as explained in Appendix B, with xo = 0. For each 
trajectory, the total time in x > 0, T+, was recorded, and the 
occupation fraction, A = T+/t, was calculated. The figure 
shows the PDF for long times of the occupation fraction A for 
a = 1,0.75,0.5,0.25 (symbols). Lamperti's PDF, Eq. (J43J) , is 
plotted as lines (For a = 1, the PDF for the simulations and 
the theory was scaled by 3 for visibility). While for a = 1, A 
is very narrowly distributed around 1/2, for a < 1, the PDF 
becomes wider and even attains a U shape for small enough 



2. An application of the occupation time functional — the 
first passage time PDF 

As a side note, we demonstrate how the fractional 
Feynman-Kac equation for the occupation time can be 
applied in an elegant manner to the problem of the first 
passage time (FPT). The FPT in the box [—§,■§] is 
defined as the time tf it takes a particle starting at 
xq = — b (0 < b < L/2) to reach x = for the first time 



41| . A relation between the occupation time functional 
of the previous subsection and the FPT was proposed by 
Kac H|: 

Pv{t f > t} = Pr{ max x(t) < 0} = lim G Xo (p,t), 

0<T<t p— >00 

where as in the previous subsection, G Xo (p, s) is the PDF 
of T + = J Q[x(t)](1t. The last equation is true since 
G Xo (p,t) = r °° e-P T +G X0 (T + ,t)dT+, and thus, if the 
particle has never crossed x = 0, we have T + = and 
e~ pT+ = 1, while otherwise, T + > and for p — > oo, 
e~ pT+ = 0. Substituting .To = — b and p — > oo in Eq. 
(UOJ) of the previous subsection gives 



lim G-b(p, s) = - < 1 

p— s-oo s 



cosh 






cosh ^ 





(45) 



the exponent is a/2, compared to a here. An equation 
for G Xo {p, s) for x$ ^ can be derived in exactly the 
same manner, leading, for long times, to Eqs. (|42[) and 
(|43]). as expected. 

For a = 1, it is easy to see from Eq. (j42]) that 
G(T + ,t) = 6(T+ - t/2) or A = 1/2. This is the expected 
result based on the ergodicity of normal diffusion. As a 
decreases below 1, the delta function spreads out to form 
a W shape. For even smaller values of a (< 0.59 [4oj|). 
the peak at A = 1/2 disappears and the PDF attains a 
U shape, indicating that the particle spends almost its 
entire time in only one of the half-boxes. For a — > 0, 
G(A) = <5(A)/2 + <5(A - l)/2, as expected. This behavior 
is demonstrated and compared to simulations in Figured] 
Details on the simulation method are given in Appendix 
B. 

For short times, (i/r) Q / 2 <C 1, we substitute in Eq. 
(|4T| the limit (sr) Q / 2 > 1, 



Go(p, s) 



s «/2-i + ( s+p )«/ 

s a/2 _|_ ( s +p)a/2 



2-1 



(44) 



In t space, this gives again the Lamperti PDF, but now 
with index a/2. This is exactly the PDF of the occu- 
pation fraction of a free particle, which is expected, be- 
cause for short times the particle does not interact with 
the boundaries Q. It can be shown that for short times, 
G X0>Q (T + ,t) = S(T + - t), and G xo<0 (T + ,t) = S(T + ), as 
expected. 



The first passage time PDF satisfies f(t) = 
[1 — Pr{£/ > <}]. We therefore have in Laplace 
space, 



/(*) 



cosh 






cosh | 





For long times, the small s limit gives 



2K n 



For < a < 1, inverting s — > t, 

fits) 



b(L - b) -(i+q) 
2K a \T(-a)\ f 



(46) 



Therefore, f(t f ) 



t 



-(!+«) 



(compared to f(tf) 



^ (i+q/2) a £ ree particle [25|, H2]), indicating that for 
a < L (tf) = oo. Eqs. (|45|) and (146)) agree with previous 
work 



3. The fluctuations 

Eq. (|41[) , giving Gq(j>, s) for the occupation time func- 
tional, cannot be directly inverted. It can nevertheless 
be used to calculate the first few moments using 



p=0 



8 



The first moment (for xo = 0) is of course (T+) 
(A) = 1/2. For the second moment, 



a(sr) 



«/2 



4s 3 2s 3 sinh [2(sr) a / 2 ] 
The long times, we take the limit of small s, 



t/2 or 



(47) 



Q 



2s 3 6s 3 ~ Q ' 

Inverting and dividing by i 2 , we obtain the fluctuations 
of the occupation fraction, ((AA) 2 ) = (A 2 ) — (A) , 

1 — a a 



<(AA) 2 ) 



6r(3 - a) 



(48) 



For a < 1 and t — > oo, we see from Eq. (|48| that 
<(AA) 2 ) = 1=2 > o. For a = 1, ((AA) 2 ) -> as t -t oo. 
The convergence to the long-times limit exhibits a t~ a 
decay. For xo ^ 0, the first moment approaches 1/2 as 
(A) ss 1/2 + 2k r^2-°a) anc ^ ^ nc fluctuations remain 
the same as in Eq. (|48| up to order t~ a . 

For short times (and xq = 0), taking the limit 
(sr) Q / 2 > 1 in Eq. g7]) gives (T 2 ) 



It-?- from which 



((AA) 2 ) 



1 - a/2 



(49) 



This is the expected result, since for short times the PDF 
is Lamperti's with index a/2 (Eq. (|4"4"|l ). 

The fluctuations ((AA) 2 ) are plotted in Figure [2] and 
agree well with Eq. (|49|) for short times and with Eq. 
(|48|l for long times. As expected, the approach to the 
asymptotic limit is slower as a becomes smaller. 



B. The time-averaged position in a harmonic 
potential 

We consider the time-averaged position, x = 
1 J* x{t)cLt, for a sub-diffusing particle in a harmonic po- 
tential, V(x) = mu> 2 x 2 /2 (fractional Ornstein-Uhlenbeck 
process [3ll.l44jh 



1. The distribution 

We first study the PDF in the long-times limit using 
the general equation (|38p . Define the second moment in 



thermal equilibrium as (^ 2 ) th = fcsT/(mw 2 ). Measuring 
x in units of */ (x 2 ) th , we have for t — > oo, 



G(x) = 



th , 



where 



s(y) 



sin(-7ra) 



£-ifo)#fo) + £-i(v)^(y) 



K(y)] 2 + K(2/)] 5 



2cos(7ra)/>( 2 /)/<(y) ; 

(50) 
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FIG. 2: The fluctuations of the occupation fraction in half 
box. CTRW trajectories were generated as explained in Ap- 
pendix B (with Xo — 0) and the occupation fraction in half 
box, A = T+ /t, was calculated. The figure shows the fluc- 
tuations ((AA) 2 ) vs. t for a = 0.4,0.7,1 (symbols). Theory 
for long times, Eq. (|48[) , is plotted as dotted lines. The fluc- 
tuations are initially equal to their free particle counterpart, 
(1 — q/2)/4 (Eq. (|49[) . indicated as dashed lines), and then 
decay to their asymptotic value, (1 — a)/4 (also indicated as 
dashed lines), as t~ a . Only for a — 1, the fluctuations vanish 
for t — > oo. 



with 



e x 2 (x - y) a dx ; 1^ 



e ^ (y — x) a dx. 



Using Mathematica, we could express the solution of the 
integrals in Eq. (|50|) in terms of Kummcr's functions. 
The full expression is given in Appendix C (Eq. (|66[) h It 
can be shown that for a = 1, G(x) = 6(x), as expected for 
an ergodic system d,[I3|. For a < 1, G(x) has a non-zero 

width, and when a — > 0, G(x) 



2k B T 

since for a — > 0, 
to (y < 1)> ff(2/) has 



27rfc fl T eX P 

which is the Boltzmann distribution 
x -> x 0, For x < 

a Taylor expansion around y = of the form g(y) 

r !jjr('$) ) + °^ 2 )- For s » (y » i), 9(v) - 

r ( a )j^ 1 ^' a ) y- 2a e~ y2 / 2 , which gives the expected results 
for a — > and a = 1. Eq. (|50p is plotted and compared 



to simulations in Fig [3] 

For s/iort times, t a -C ( a;2 ) th /-Kaj the particle is at the 
minimum of the potential and therefore behaves as a free 
particle. For the free-particle case, we have previously 
shown a scaling form for Xq = [25j 



G(x, t) 



1 



fKlt a / 2 



fKZt a / 2 



(51) 



where h a (y) is a dimcnsionlcss scaling function. This 
behavior is numerically demonstrated in Fig [3] 
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FIG. 3: The PDF G(x) for a particle in a binding harmonic 
field. CTRW trajectories were generated using the method 
described in Appendix B, with xo = 0. Top panel: Simula- 
tion results for long times for a = 0.25, 0.5, 0.75, 1 (symbols). 
Theory for t — > oo, Eq. (|50[) . is plotted as solid lines (For 
a = 1, the PDF for the simulations and the theory was scaled 
by 2 for visibility). For a = 1, the distribution is a delta func- 
tion, whereas for a < 1, x is a random variable even for long 
times, indicating ergodicity breaking. Bottom panel: Simu- 
lation results for the PDF of x for a number of short times 
and for a = 0.25,0.5,1 (symbols). The plot illustrates the 
free-particle scaling form, Eq. (f5lT) . 



2. The fluctuations 

The PDF of the time-averaged position was shown in 
the previous subsection to have a non-trivial limiting 
distributions for t -> oo (Eq. ((50)) ) and t — S- (Eq. 
(f5Tj) ). However, the shape of the PDF for other times 
is unknown. In this subsection, we show that using the 
fractional Feynman-Kac equation, we can determine the 
width of the distribution for all times. 

Let us write the forward equation in (p, s) space for 
the functional A = xt = f* x(r)dr (U(x) = x) and for 
xq = 0. Since A is not necessarily positive, p here is the 



Fourier pair of A and we use Eq 
sG(x,p, s 



if- 



5{x) = ipxG(x,p, s 
d mw 2 ! 1 



of Section [EEO 
(52) 



+ 



dx 2 dx k B T 



s — ipx] 1 a G(x,p,s). 



To find {A 2 ), we use the relation 



f°° B 2 



dx. 



p=0 



Operating on both sides of Eq. ([52]) with — -§^, substi- 
tuting p = 0, and integrating over all x, we obtain, in s 
space, 



;(A 2 ) =2{Ax) i 



(53) 



where we used the fact that the integral over the Fokkcr- 
Planck operator vanishes. Eq. (|53[) can be intuitively 

understood by noting that {A 2 ) = 2 ^AA^ and A = x. 

We next use Eq. ([52]) and 



f°° d 

( Ax ) s = ~ l x ^r G (x,p,s) 
, op 



dx, 



p=0 



to obtain, 

s (Ax), = [1 + (1 - a)(sr)- a ] (x 2 ) s - s(st)-« (Ax), , 

where we defined the relaxation time r Q = 
k B T/(K a muj 2 ) = (x 2 ) th /K a . Thus, 

Finally, to find (x 2 ) s , we use (x 2 ) s = J^, x 2 G(x,p = 
0, s)dx, 

S (x 2 ) s = 2X qS - q -2 S ( S t)- q <x 2 ) s , 

where we used the normalization condition J G(x,p — 
0,s)dx = 1/s. Thus, 



2<- 2 ) th 



(55) 



2 + (st) q ' 

Combining Eqs. ([55]). ([5I|). and ([55]). we find, 

_ 4 (1 - q) + (sr)° (^ 2 ) th 
^ ' s s 3 l + (sr) Q 2 + (sr) Q ' 

To invert to the time domain, we write {A 2 } g as partial 
fractions: 



(l-a) + 2a 7 i^-(l + «) ; (ST)Q 



(56) 



l + (sr) Q y 2 + (sr) Q 

Inverting the last equation, we find 

(A 2 ) = (x 2 ) th t 2 x (57) 
{(1 - a) + 4aE a , 3 [~(t/r) a ] - 2(1 + a)E a , 3 [-2(t/r) a ]} , 
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where we used the Laplace transform relation [45| 

1 (st)° 



r 

Jo 



e- st t A E a , 3 [-c{t/T) a ]dt 



i s 3 c+{sr) a ' 
and E a3 (z) is the Mittag-Lcmer function, defined as [45 



E a , 3 (z) = 



r(3 + an) 



To obtain the fluctuations of the time-averaged position, 
({Ax) 2 ) = (x 2 ) - (x) 2 , we use (x 2 ) = (A 2 ) /t 2 and 
(x) = (since Xq = 0). This gives 

((Ax) 2 ) = (x 2 ) th x (58) 
{(1 - a) + AaE a , 3 [~(t/r) Q ] - 2(1 + a)E a , 3 [-2(t/r) a }} . 



Eq. (|58|1 is plotted (using [46j|) and compared to simula- 
tions in the top panel of Figure 01 

To find the long times behavior of the fluctuations (1551) , 
we expand Eq. (|56p for small s, invert, and divide by t 2 , 



(3a-l)(x 2 \ h ft 
T(3-a) 



((Ax) 2 ) « (l-a){x 2 ) + 



(59) 

Thus, for a < 1 and i -> oo, ((Ax) 2 ) = (1-a) (z 2 ) th > 
and ergodicity is broken. Only when a = 1, we have 
crgodic behavior ((Ax) 2 ) = 0. As we observed for the 
occupation fraction (Eq. (|48|) ). Eq. (|59|) too exhibits a 
t -Q convergence of the fluctuations to their asymptotic 
limit. 

For short times, 



E, 



Q,3 



(t/r) a 
r(3 + a)' 



Therefore, 



((Ax) 2 ) 



4(x 2 ) 



th 



r(3 + a) 



(60) 



Noting that (x 2 ) th /T a = K ai we can rewrite Eq. (|60|). 
as ((Ax) 2 ) w r ^° a ^ °, which is, as expected, equal to 

the free particle expression (25[. 

The bottom panel of Figure 2] presents the fluctua- 
tions of the time-average (for xq = 0) for a wide range 
of times and for a = 0.05, 0.1, 0.15, 1. As expected 
from Eqs. (|5"9"|) and (|60[) . the fluctuations increase from 
((Ax) 2 ) = at t —> to their asymptotic value at t —> oo, 
(x 2 ) th (1 — a). However, as can be seen also in Eq. (|59|) . 
for a > 1/3 the fluctuations display a maximum and 
decay to their asymptotic limit from above. We found 
numerically that the value of the maximal fluctuations 
scales roughly as a -1 / 2 (not shown). It can also be seen 
that for almost all times and all values of a, the fluctu- 
ations ((Ax) 2 ) decrease as the diffusion becomes more 
"normal" (increasing a), as expected. However, this pat- 
tern surprisingly breaks down for a < 0.15, for which 




FIG. 4: The fluctuations ((Ax) 2 ) for a particle in a har- 
monic potential. Top panel: CTRW trajectories were gener- 
ated using the method described in Appendix B, with xo = 0. 
Symbols represent simulation results for a — 0.25, 0.5, 0.75, 1. 
Theory, Eq. (|58jl. is plotted as solid lines. The straight 
dashed lines are limt-,.^ ((Ax) 2 ) = (1 — a) (x 2 ) th . Except 
for a = 1, the fluctuations do not vanish when t — >• oo and 
thus ergodicity is broken. The dotted lines represent the 
long-times and short-times approximations, Eqs. ()59|) and 
(|60[) . respectively. Bottom panel: The fluctuations, Eq. (|58[) . 
plotted for a wide time range [10 -15 , 10 30 ]. Shown are 20 
curves for a = 0.05, 0.1, 0.15, 1 (top to bottom). The fluc- 
tuations display a maximum when a < 1/3 and a crossover 
when a < 0.15. As expected, the fluctuations approach their 
asymptotic value slower for smaller values of a. 



there is a time window when the fluctuations increase 
with a. 

It is straightforward to generalize our results to any 
initial condition with first moment (xo) and second mo- 
ment (xq). The first moment of the time-average is 
(x) = (xq) E at 2 [— (t/ r ) a ]i which decays for long-times 
as (x) ~ r (2-a) (t) ^ e secon d moment is 
(x 2 ) = (1 - a) (x 2 ) th + 2a [2 (x 2 ) th - (x 2 )] E a . z [-{t/rT} 
+ 2(1 + a) [(x 2 ) - (x 2 ) t J £ Q , 3 [-2(t/r) a ] , (61) 

from which the fluctuations directly follow. For long 
times, 

((Ax) 2 )«(l-a)(x 2 ) th 

(3a-l)(x 2 ) th + (l-a)(x 2 ) ft 



T(3 - a) 
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For short times, 



<(Ax) 2 )«((Ax ) 2 )-2 



<(Az ) 2 ) 2<x 2 ) 



th 



T(2 + a) T(3 + a) 



where ((Axo) 2 ) = (x 2 ) ) — (xo) 2 . According to the last two 
equations, if the system is already in equilibrium at t = 
such that (xq) = (x 2 ) th , the fluctuations monotonically 
decay, for all a, from (s 2 ) th at t = to (x 2 ) th (1 — a) at 
t — > oo. 

For a = 1 (and = 0), we find the known result 47j: 



((^) 2 ) Q= i = 7 



4e 



-t/7 



s -2t/r + 2* _ g 



(62) 

To derive the last equation, we used the relation 
£■1,3(2) = [e z — z — l]/z 2 - Since the ordinary (a = 1) 
Ornstein-Uhlenbeck process is a Gaussian process 
the PDF of 2; is a Gaussian too, with the variance indi- 
cated by Eq. ([62)) . 



3. Fractional Kramers equation 

Finally, we remark on the connection between the frac- 
tional Feynman-Kac equation of this subsection and an 
important class of processes in which the velocity of the 
particle is the quantity undergoing sub-diffusion. For 
such processes, Friedrich and coworkers have recently 
developed a fractional Kramers equation for the joint 
position- velocity PDF (2(| |27|- For example, consider 
a Rayleigh-like model in which a free, heavy test particle 
of mass M collides with light bath particles at random 
times, but where the times between collisions are dis- 
tributed according to ip(r) ~ t"( 1+q ). The PDF of the 
velocity of the test particle, GJv,t), satisfies the frac- 
tional Fokker-Planck equation [44[ : 



d_ 

dt 



G(«,t) = 7o 



k B T d 2 



d 



M dv 2 dv 



where Drl" ls the Riemann-Liouville fractional deriva- 
tive operator (sec Section III A3[) and 7 Q is the damping 
coefficient. Since in the collisions model x (t) = J v(t)cLt, 
x is a functional of the trajectory v(t), and therefore, 
the joint PDF of x and v, G(v,x,i), is described by our 
fractional Feynman-Kac equation. Denoting the Fourier 
transform x — > p of G(v,x,t) as G(v,p,t), we have (see 
Eq. (HID), 



d_ 

dt 



G(v,p,t) 



+ 7 



ipvG(v,p, t) 
k B T d 2 



d 



M dv 2 dv 



(63) 
V\-<*G{v,p,t), 



model, for < a < 1 the motion is ballistic, (a; 2 ) ~ t 2 , 
while for a = 1 it is diffusive, (x 2 ) ~ t (see Eq. (f59|0 . Eq. 
(|63p is exactly equal to the fractional Kramers equation 
derived by Friedrich and coworkers [2^, [27j , and in that 
sense, our results generalize their pioneering work. 



V. SUMMARY AND DISCUSSION 

Time-averages of sub-diffusive continuous-time ran- 
dom walks (CTRW) in binding fields are known to exhibit 
weak ergodicity breaking and were thus the subject of re- 
cent interest. In this paper, we used the Feynman-Kac 
approach to develop a general equation for time-averages 
of CTRW (Eq. (|22p). which can be seen as a fractional 
generalization of the Feynman-Kac equation for Brown- 
ian motion. The equation we derived describes the dis- 
tribution of time-averages for all observables, potentials, 
and times. We also derived a backward equation (Eq. 
(|3ip) which is useful in practical problems. 

We investigated two applications of our equations: the 
occupation fraction in the positive half of a box, and 
the time-averaged position in a harmonic potential. In 
both cases, we obtained expressions for the PDF for long 
times and for short times and calculated the fluctuations. 
We found that the fluctuations decay as t~ a to their 
asymptotic limit, which is non-zero for anomalous diffu- 
sion, a < 1. Our fractional Feynman-Kac equation thus 
provides a general tool for the treatment of time-averages 
and for the study of the kinetics of weak ergodicity break- 
ing. 

Recently, the occupation time functional has been 
studied in the context of dynamical systems with an in- 
finite (non-normalizable) invariant measure (49j . It re- 
mains to be seen whether a framework similar to that 
of the fractional Feynman-Kac equation could be devel- 
oped for general functionals of these processes. We also 
note that while the (integer) Feynman-Kac equation can 
be derived using path integrals (22[, a path integral ap- 
proach for functionals of anomalous sub-diffusion is still 
awaiting (but see preliminary results in the upcoming 
book |50j). 
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Appendix A: Time-dependent forces 



where D t a is the fractional substantial derivative, here 
equal in Laplace s space to (s — ipv) 1 ^". Within this 



In our model of CTRW with a time-dependent force, 
jump probabilities are determined according to the force 
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at the time of the jump. To derive an equation for 
G(x, A, t) in that case, we rewrite Eq. ([7]) as follows: 



X (x,A,t) = G (x)S(A)S(t) 



(64) 



+ / ip(r)L(x + a,t)x[x + a, A - tU(x + a),t - T]dr 



+ / i/)(t)R(x — a, t)x[x — a, A — tU(x — a), t — rjdr. 
Jo 

Note that the jump probabilities are time-dependent (but 
have no memory) . Laplace transforming A — » p and t — > 
s, using the Laplace identity L{tf(t)} = — Jj/(s), 

X(x,p,s) = G Q (x) 

( d \ - 

+ L ix + a, ~-q^J i-A s +pU(x + a)]x{x + a,p 7 s) 



R[x-a, -— ) ijj[s+pU(x - a)]x(x - a,p,s). 



Fourier transforming x — > k, 



X(k,p,s) = G (k)- 



cos(fca) + i sin(fca 

d_ 

dk 



s + pU —i 



2k B T 
X(k,p, s). 



Continuing as in Section HTA21 we find the formal solu- 
tions for x(fc,p, s) and G(k,p,s) and then take the con- 
tinuum limit. This gives: 

sG(k,p, s) ~ Go(fc) = -pU (-i-^A G(k,p, s) 



ik 



S + pU \- l dk 



k B T 



G(k,p, s) 



In this limit, a — >• and B a — >• but the generalized 
diffusion coefficient K a = a 2 /{2B a ) (Eq. (|2T])) is kept 
finite [28]. We simulate trajectories of this process as 
follows [51j. We place a particle on a one-dimensional 
lattice in initial position xq, where usually xq = 0. We 
set the lattice spacing a and the generalized diffusion 
coefficient K a and determine B a = a 2 /(2K a ). Wait- 
ing times are then drawn for a = 1 from an exponen- 
tial distribution ^>{t) = e~ T / T °/ro with mean tq = B\. 
This is implemented by setting r = — roln(u), where u 
is a number uniformly distributed in [0,1]. For a < 1, 
we set t = [B a /T(l — a)] 1 /" and r = t u 
corresponds to ip{r) = | r ^_° a ^| T~t 1+ ") (r > 
([5])). After waiting time r, we move the particle right 
or left with probabilities R(x) or L(x), respectively, as 
given by Eq. (HU). For the harmonic potential, Eq. ([14")) 

gives R(x) = ^ (l 
Since the typical 



~ 1 / Q , which 
tq; see Eq. 



2(x 2 



andL(.x) = | (l 



2(.i 



is of the order of y / [x 2 



it is 



sufficient to choose a <C \/ (x 2 ) th to guarantee that 
< R(x),L(x) < 1 (sec discussion in [43). For the 
box, R(x) = L(x) = 1/2 and we make the boundaries at 
x = ±-j reflecting. 

The parameters we used in the simulations were as fol- 
lows. In all simulations, wc used a = 0.1 or smaller, and 
each curve represents at least 10 4 trajectories. For the 
occupation time in a box, we set L = 2 and K a = 1 , and 
the final simulation time in Figure [1] was t = 10 3 . For 
the time-averaged position in the harmonic potential, we 
set K a = 1/2 and (x 2 ) th = 1/2 (or r Q = 1). In Fig- 
ure [U the final simulation times were as follows. For the 
long-times limit (top panel) we used t = 10 7 , 10 4 , 10 3 , 10 3 
for a = 0.25,0.5,0.75,1, respectively. For the short 
times (bottom panel), we used t = 10 -3 , 10 -2 , 10 _1 
for a = 1, t = 10~ 5 , 10~ 4 , 10 -3 for a = 0.5, and 
t = 10" 6 , 10- 5 , 10~ 4 for a = 0.25. 



Inverting k — > x, s — > t, we obtain the fractional 
Feynman-Kac equation for a time-dependent force: 

^G{x,p,t) = K a L$lvl- a G( X ,p,t) - P U(x)G(x,p,t), 

(65) 

where 



(t ) _ d 2 8 F(x,t) 



FP 



dx 2 dx k B T 
is the time-dependent Fokkcr-Planck operator. 



Appendix C: The t — > oo distribution of the 
time-averaged position in a harmonic potential. 



Consider the time-averaged position, x = | J * x(t)cIt, 
for a sub-diffusing particle in a harmonic potential, 
V(x) = mu! 2 x 2 /2. Using the thermal second moment, 



• 2 ) - 

/th 



k B T I '(raw 2 ), and for t — > oo, we have 



Appendix B: The simulation method 



G{x) 



'th, 



The fractional Feynman-Kac equation describes the 
joint PDF of x and A in the continuum limit of CTRW. where 
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. . sm(7ra) 

g(y) = x 



y r - r(i + q) 



y 2 \ v ( 1 , a 3 y 2 



-V2T(a)T 



)2+a 2t^2 



1 + a 



y^ l 



e^M 2 



^V^e^VTTJ/rCl + a)M 
'1 + a 



2 ' 2' 2 
1 + a 3 y 2 
2 ' 2' Y 
1 - a 3 y 2 
2 '2'~"2~ 
1 - a 3 y 2 



U 1 



2' 2' 2 
a 3 if 
2' 2' 2 



-2M 
2M 



1-a 3 y 2 
2 '2'~~2 
1 + a 1 y 2 



U 



- 2cos(7ra)M 2 1 + 



2 ' 2' 2 

a 3 y^ 
2' 2' 2 



U 



a 1 y_ 2x 
2' 2' 2 

a 1 y 2 / 
2' 2' 2 



w M 2 



2 ' 2' 2 
2' 2' 2 



M 



a 1 jp_ 

'2'2' 2 



2cos(7ra)M 2 



1 + a 1 j/ 



2 \ 1 



W(1 + o)D . |1 + «|,£ 



(66) 



In the last equation, M(et, b, z) is the confluent hyper- 
geometric (or Kummer's) function of the first kind and 
U(a, b, z) is the confluent hypergeometric (or Kummer's) 



function of the second kind (52|. Eq. (|6"6"|) is valid 
for y > 0. Due to the symmetry of the potential, 

g(-y) = g{y)- 
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